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Abstract 

DNA torsion dynamics is essential in the transcription process; simple 
models for it have been proposed by several authors, in particular Yaku- 
shevich (Y model). These are strongly related to models of DNA sepa- 
ration dynamics such as the one first proposed by Peyrard and Bishop 
(and developed by Dauxois, Barbi, Cocco and Monasson among others), 
but support topological solitons. We recently developed a "composite" 
version of the Y model, in which the sugar-phosphate group and the base 
are described by separate degrees of freedom. This at the same time fits 
experimental data better than the simple Y model, and shows dynamical 
phenomena, which are of interest beyond DNA dynamics. Of particular 
relevance are the mechanism for selecting the speed of solitons by tuning 
the physical parameters of the non linear medium and the hierarchal sep- 
aration of the relevant degrees of freedom in "master" and "slave" . These 
mechanisms apply not only do DNA, but also to more general macro- 
molecules, as we show concretely by considering polyethylene. 
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Introduction 



Following the early works of Davydov on solitons in biological systems [13], 
it has been conjectured since a long time [17] that nonlinear excitations - in 
particular, kink solitons or breathers - could be present in the DNA double 
chain and could play a functional role, in particular in the processes of DNA 
denaturation and transcription. 

This general idea meets of course essential difficulties when one tries to 
translate it into quantitative terms due to the formidable complexity of the DNA 
molecule [9, 45]. This is organized in two helices; each of them is composed of 
adjoining nucleotides. A nucleotide consists of a unit of the sugar-phosphate 
backbone (identical in each nucleotide) and an attached nitrogen base (this 
can be of four different types; the base at a given site on one chain uniquely 
determines the base at the same site on the other chain, as they must be one of 
the Watson-Crick pairs). This makes 30-35 atoms and hence about 100 classical 
degrees of freedom for each nucleotide; each helix is then made, depending on 
the species, of 10 6 — 10 9 nucleotides. 

In view of the quasi-regular structure of DNA - and despite the fact genetic 
information is embodied in the non-regular part of the structure - it is quite 
reasonable to start the modelling by considering a polymer made of identical 
units (the nucleotides), deferring taking into account the actual base sequence 
and hence inhomogeneities in the structure to a later moment (and to computer 
simulations rather than analytical investigation). 

Needless to say, DNA like any other molecule actually obeys quantum rather 
than classical mechanics. The first consequence of this is that nucleotides can be 
realistically thought as made of rather rigid subunits, and one can just consider 
the degrees of freedom of these subunits [9, 26, 45]. Under closer scrutiny, it 
turns out that some of these degrees of freedom are more easily excited and 
hence dominant, i.e. those related to a radial movement of the bases away 
from the double helix axis, and those related to rotations of the bases and the 
sugar ring in a plane nearly orthogonal to the double helix axis; in the standard 
nomenclature of DNA deformations [9, 45, 26], they correspond respectively to 
stretch and opening. 

These considerations are at the basis of DNA modelling as considered in 
the Nonlinear Mathematics and Theoretical Physics communities, where one 
aims at reproducing significant experimental observations on the basis of models 
with few degrees of freedom per nucleotide. These cannot substitute for more 
massive quantum chemistry computations, but could identify relevant degrees 
of freedom, hence help in organizing our understanding of the complex DNA 
dynamics. It is worth stressing, in this respect, that DNA is not only complex 
structurally, but also performs a great wealth of biological tasks. It is thus 
not impossible that one can consider different models of it depending on the 
biological process one aims at modelling; from this perspective, considering only 
a few degrees of freedom per nucleotide in a model aiming at a specific mode 
of DNA dynamics, relevant in a specific process, is quite reasonable despite the 
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Figure 1: A schematic view of the transcription process. RNA Polymerase reads 
the bases sequence and produces RNA messenger; the reading requires local 
unwinding of DNA double helix in a "transcription bubble" region involving 
about 20 base sites. When active in transcription, RNA Polymerase travels 
along DNA at a speed of about 0.5 — 1.0 x 10 3 bases per second, and the 
transcription bubble moves along the molecule at the same speed. Adapted 
with modifications from [9]. 

underlying overall complexity of the molecule and its spatial organization. 1 

As mentioned above, simple DNA models are primarily focusing on DNA 
denaturation and transcription; more specifically, they aim at describing the 
deformation of DNA structure corresponding to the two main degrees of freedom 
mentioned above, "radial" ones (movement of the bases directly away from the 
double helix axis), thought to be relevant in DNA denaturation; and "torsional" 
ones (rotations of the bases and the sugar ring in a plane orthogonal to the 
double helix axis), thought to be relevant in DNA transcription where local 
untwisting of the double helix would make possible to RNA Polymerase to 
access the base sequence without disrupting the double helix [9, 17, 26, 45, 54]; 
see Figure 1. 

In this note we will first discuss these models, in particular some concrete 
models widely studied in the literature, with their success and limitations; and 
then consider a specific concrete "composite" model we recently proposed for 
the torsional dynamics and which on the one hand is free of limitations pertain- 
ing to other torsional models, and on the other hand shows phenomena which 
are obviously not specific to it, and could be of interest - theoretical but also 
applicative - in a much wider nonlinear mechanics context. In order to illustrate 
this wider applicability, in section 5 we also discuss how these phenomena show 
up in a different macromolecule, i.e. polyethylene. 2 

lr The models we will consider take into account only the double helical structure, disregard- 
ing the way this double helix is organized in three-dimensional space; that is, we are actually 
focusing at the DNA structure on small length scales. 

2 Our discussion of the DNA composite model will go over our recent results [5, 6, 7, 8]; on 
the other hand, the discussion of how the general mechanisms devised in that context apply 
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We will not discuss DNA structure and its Physics, for which the reader is 
referred e.g. to [9, 45] and [19, 58] respectively. Experimental studies of single- 
molecule DNA dynamics are discussed e.g. in [35, 43, 48]; for general discussion 
of modelling DNA at different scales, see [4]. 

1 Mechanical models of DNA 

A number of mechanical models of the DNA double chain have been proposed 
over the years, focusing on different aspects of the DNA molecule and on dif- 
ferent biological, physical and chemical processes in which DNA is involved. A 
discussion of such attempts is given in the books by Yakushevich [54] and by 
Dauxois and Peyrard [42] , as well as in the review paper by Peyrard [39] . 

In recent years, two models have been extensively studied in the Nonlin- 
ear Physics literature. These are the "radial" model by Peyrard and Bishop 
[40] (and the extensions of this formulated by Dauxois [12], Dauxois, Peyrard 
and Bishop [41], and later on by Barbi, Cocco, Peyrard and Ruffo [1, 2]; see 
also Cocco and Monasson [11]. More recent advances are discussed in [39] and 
references therein); and the "torsional" one by Yakushevich [51], which had pre- 
cursors discussed in [54] and is put in perspective within a hierarchy of DNA 
models in [52, 53]. We will refer to these as the PB and the Y models respec- 
tively. 3 

The interplay between radial and torsional degrees of freedom of bases is 
considered organically in the Barbi-Cocco-Peyrard (BCP) model [1, 2, 11]; the 
latter was formulated as an extension of the Peyrard-Bishop-Dauxois model 
[40, 41], i.e. in the context of "radial" dynamics. 

1.1 The PB model 

In PB-likc models, the bases can only move radially away from the double helix 
axis. The potential energy corresponds to stacking interactions between succes- 
sive bases on each chain, and pairing interactions between bases at corresponding 
sites on opposite chains [9, 26, 42, 45]. 

We denote by G 1Z + the position of the base at site n on the helix (±). 
It is convenient, for our present purposes, to limit discussion of this and other 
models to symmetric configurations, i.e. configurations withri +) (t) =r [ n \t) = 
r n (t) at all times. In this case, the Lagrangian describing the PB model is (with 

to polyethylene — based on the Zhang-Collins model [62] - is new. 

3 Here we will discuss, for the sake of brevity, only the "planar" versions of these models, i.e. 
overlook so called "helicoidal" interactions. These are interactions between bases which are not 
first-neighboring in the chain sequence but which come to be near in three-dimensional space 
due to the helical geometry of the DNA molecule. Considering these introduce qualitative 
differences in the dispersion relations, both in the PB [12] and in the Y model [20]; sec [22, 26] 
for a discussion. The same will apply for the composite Y model discussed in Sect. 2; see [5] 
for its full version, taking into account "helicoidal" interactions. 
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m the mass of bases) 



£pb = ^2m(r n ) 2 - K s (r n+1 -r n f - V p (r n ) . (1) 

n 

In the PB model, V p is a Morse potential with a minimum at the equilibrium 
position, corresponding to r n = p ~ 2A: 

K p (r) - D (cxp [-a(r - p)] - l) 2 . (2) 

The Euler-Lagrange (EL) equations for £pb are of course 

mr n = if s (r„ + i - 2r„ + r„_i) - ^V p '(r n ). (3) 

We will pass to consider the continuum approximation for these, i.e. substitute 
the infinite array of scalar variables {r n (t)} with the interpolating field R(x, t), 
such that R(nS,t) w r„(i); here (5 ~ 3.4A is the distance between successive 
base pairs positions along the double helix axis. Using now a second order 
approximation, 

R(x ± J, t) M R(x, t) ± ^(x, i) + (5 2 /2)R xx (x, t) , (4) 

and writing for short k = K S S 2 , the (3) yield the nonlinear wave equation 4 

mRtt = kR xx - (l/2)[dV(R)/dR\ . (5) 

It can be shown that this supports breathers; their size and oscillation fre- 
quency (choosing parameters so that V p describes as far as possible a Hydrogen- 
bond interaction) are compatible with those observed in real DNA. The discrete 
model can be put in a thermal bath and numerically simulated, and again one 
observes a behavior compatible with the one observed in DNA denaturation 
[39, 42]. More refined versions of this model [41] consider modified expressions 
for the stacking energy, improving - also qualitatively - the correspondence with 
actual DNA behavior; we refer again to [39, 42] for details. 

1.2 The Y model 

In the simple models for DNA torsional dynamics, one studies a system of non- 
linear equations which in the continuum limit reduce to a pair of sine-Gordon 
(SG) type equations; the relevant nonlinear excitations are kink solitons - which 
are solitons in both dynamical and topological sense - which describe the un- 
winding of the double helix in a "bubble" . 

The main biological interest of these model lies in the identification of this 
unwound bubble with the transcription region (this is indeed an "open bubble" 

4 The same equation is also obtained passing to the continuum approximation directly in 
the Lagrangian, i.e. considering the Lagrangian density Lpb = Rt ~ KsR% — V P (R). The 
same holds for the other models considered below. 
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of about 20 bases, to which RNA Polymerase (RNAP) binds; the RNAP travels 
along the DNA double chain, and so does the unwound region). The idea of 
Englander et al. [17] was that the open bubble could correspond to nonlin- 
ear excitations and thus be present due to the nonlinear dynamics of the DNA 
double helix itself; the RNAP would then use them to travel along DNA. In 
this way a number of questions - in particular, concerning energy flows - would 
receive a simple explanation. Note that their model, and subsequent ones con- 
tinuing their research, are not concerned with the DNA-RNAP complex, but 
the dynamics of the DNA double helix alone. 

To be specific, let us consider the model proposed by Yakushevich (see e.g. 
[54] for similar models proposed earlier on by other authors, starting with [18, 
50, 56]); now the degrees of freedom for the rotation of the base at site n on 
the chain (±) will be denoted as ipit\ and again we restrict to the symmetric 
case, i.e. enforce ip\t\t) = tp^~\t) = ip n (t) at all times (this makes that in 
the continuum approximation we will get a single equation of SG type rather 
than two coupled ones). The Lagrangian describing the Y model is (here m is 
a moment of inertia) 

£y = ^2m(ip n ) 2 - K s - %(<Pn) 5 (6) 

n 

the choice by Yakushevich for the intrapair potential was that of a simple har- 
monic potential 5 [51], resulting in (we stress here r is a geometrical constant, 
not a dynamical variable!) 

%{ip) = —AK p r 2 cos(^) . (7) 
The EL equations for Cy are a set of sine-Gordon coupled equation, 

rrvfin = K s (p n+1 -2<f n + Lp n - 1 ) - (1/2) Vp(<p n ) . (8) 

Passing to continuum approximation with interpolating field &(x,t) (where 
$>n5, t) w <p n (i), like above) and second order approximation 

$(z ± S, t) w $(x, t) ± 5$ x (x, t) + (5 2 /2)$ xx (x, t) , (9) 

the (8) reduce to a sine-Gordon equation 

m$« = k<S> xx - A sin($) ; (10) 

here of course we have made use of the explicit form of V p , set A = 2K p r 2 , and 
defined k as above. 

5 The Y model also sets to zero the equilibrium length for this harmonic interaction (con- 
tact approximation); this results, as observed by Gonzalez and Martin-Landrove [31], in a 
degeneration of the model (it is thanks to this that we obtain exactly the SG equation). If we 
go beyond the contact approximation, the equations we obtain are more complex, dispersion 
relations change quantitatively and qualitatively, but soliton solutions are little affected by 
this [23]; see also Figure 4 below. 
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As well known, the sine-Gordon equation supports topological soliton so- 
lutions [16]; these solutions are also solitons in the dynamical sense [10]. The 
basic soliton solution with speed v is 



$(ar,t) 



4 arctan [exp(/3(x — vt))] , 



(11) 



where we have written 



= 



2K v r 2 



jj, = m 2 v 2 — K S S 2 . 



(12) 



Note that v is a free parameter, subject to the condition 6 \v\ < v* = \ / K s (S/m). 
A selection of the speed based on energy is also not present, as the soliton energy 
has a very weak dependence on its speed except for v ~ ±t>„, see [21]. 

It has been shown that the Y model gives a correct prediction of quantities 
related to small amplitude dynamics, such as the frequency of small torsional 
oscillations; and also of quantities related to fully nonlinear dynamics, such as 
the size of solitonic excitations describing transcription bubbles [26, 54]. 

On the other hand, the Y model is not capable of providing a satisfactory 
prediction for other quantities: in particular, if we try to fit the observed speed 
of transversal waves along the chain [54], this is possible only upon assuming 
unphysical values for the coupling constants [55]. That is, with a physical value 
of the parameters - in particular, for K s w 120KJ/mol - the Y model predicts 7 
a speed v s=s 320m/s, while in order to get a speed of the order of 2 Km/s as 
observed in experiments 8 and [32, 55] one has to take K„ w 6000. 

Finally, we stress once again that the Y model, like the PB one and unlike 
the BCP one, assumes that there is a single (angular in this case) degree of 
freedom for each nucleotide. 



In the BCP model [1, 2, 39, 42], the state of each base is described by both 



6 The existence of a limit speed for travelling waves (and not just soliton solutions) is related 
to the Lorenz invariance of the SG equation; see also [7] and Sect. 3 below. 



7 It follows from the discussion in [55] that v = A\/K s with A (» 28.3335; see in particular 
the non- numbered formula before eq.(7) there. 

8 More precisely, measures on DNA fibers in the B-DNA conformation give v = 1.9Km/s 
[32], while measures in DNA crystals yield a speed of 2.45 Km/s, which can grow up to 4.15 
Km/s depending on counterions concentration and chemical nature [59]. Note only transverse 
waves speed matters here, as the model does not allow for longitudinal waves. 



1.3 The BCP model 
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with V p the same Morse potential as above (in the non symmetric case, pairing 
would depend on <j> variables as well) and L the equilibrium distance between 
bases in three-dimensional space 9 . As for the last term in (13), this is a curvature 
term whose role is to avoid "zig-zag" configurations, possible for Go = 0; it 
should thus be omitted in the continuum approximation. 

In this case, proceeding as above (or more simply using the continuum ver- 
sion of the Lagrangian) one obtains [28, 39, 40] in the continuum approximation 
the EL equations 

mRu = mm\ - 7 (i?p2 -R xx )- (1/2)V'(R) , , . 

mi? 2 $ tt = -2mRR t $ t + j(R 2 <t> xx + 2RR X <S> X ) . 

It has been shown that for L > S this supports nontrivial topological soliton 
solutions; the determination of these reduces to the study of a single equation 
thanks to an exact conservation law (the model is invariant by global rotations), 
but solitons are only determined numerically [28]. 

The BCP model is quite successful in describing breather excitations of the 
DNA double chain, and its predictions (or those of the cognate CM model) fit 
well experimental observations [39, 40] related to the DNA denaturation process. 

2 Composite model of DNA torsional dynamics 

A closer look at nucleotides structure and conformation [9, 45] shows that tor- 
sional motions can take place both as a rotation of the nitrogen base with respect 
to the sugar ring and as a rotation of the sugar-phosphate group 10 ; it should 
be stressed that the former one is subject to steric hindrances, i.e. is limited by 
interactions with other parts of the DNA molecule. 

Thus, one should consider "composite" torsional models, in which we de- 
scribe the state of each nucleotide by two independent angular degrees of free- 
dom, one related to the sugar-phosphate group and one to the nitrogen base; see 
Figure 2 for details. Note that in this scheme one is not considering "radial" 
(stretch) motions. 

2.1 The model 

Such a model, described graphically in Figure 2, was recently put forward and 
studied by three of us [5] (see also [6] for a discussion of it focusing on its 
mathematical features); it shows some phenomena - to be discussed in later 
sections - which are not specific to it but apply to a more general class of 
models and could be of interest in fields quite far from DNA dynamics [7, 8]. 

9 The geometry of the BCP model assumes the distance h n between successive base pairs 
planes (at sites n and n + 1) is constant and equal to S, while the length £„ of the sugar- 
phosphate backbone unit connecting them can vary. A similar model where l n is fixed and 
h n can vary has been formulated by Cocco and Monasson (CM model) [11]. 

10 More precisely, of the sugar ring around the P — O — C — C — C — O — P... chain. 
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Figure 2: A base pair in the composite Y model (reproduced from [5]). The ori- 
gin of the coordinate system is in O; the angles 9\ and #2 correspond to torsion 
of sugar-phosphate backbone with respect to the equilibrium B-DNA confor- 
mation; the angles ip\ and ip2 between A'B' and B'C correspond to rotation of 
bases around the C — N bond linking them to the nucleotide. All angles arc in 
counterclockwise direction. 

Our model can be considered as an extension of the Y model. It turns out 
that the simple Y model captures to a large extent the essential features of 
the nonlinear dynamics of the composite model. On the other hand, the more 
realistic geometry of the composite model yields a relevant improvement of the 
descriptive power of the model at both the conceptual and the phenomenological 
level; in particular, the composite Y model allows for a more realistic choice of 
the physical parameters. 

The different degrees of freedom we use will play a fundamentally different 
role in the description of DNA nonlinear dynamics. The backbone degrees of 
freedom (recall the rotations they describe are not limited) are "topological" 
and play to some extent a more relevant role, in that the solitons are mainly 
associated to them; while those associated to the base (recall the associated 
rotations are subject to steric hindrances) are "non topological" and represent 
small oscillations. These different roles are specially clear when we consider the 
limit in which our model reduces to the standard Y model, in which only the 
topological degrees of freedom are present. 

It should be stressed that this feature is specially interesting in connection 
with the possibility (discussed elsewhere [14]) to consider more realistic models, 
in which differences among bases are properly considered, as perturbations of 
our idealized uniform model. As the essential features of the fully nonlinear 
dynamics are related only to backbone degrees of freedom, such a perturbation 
should be expected to show the same kind of nonlinear dynamics as our uniform 
model. 
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The Lagrangian defining the composite Y model will be written as 



C cY = T - (U t + U s + U p ) - U c 



(15) 



with T the kinetic energy, U t is the backbone torsional potential, U s the stacking 
potential, and U p the pairing potential. 

Moreover U c is a constraining potential which represents the steric hin- 
drances to the base rotations (not accompanied by sugar-phosphate rotations); 
its explicit expression is to a large extent arbitrary, provided <p n are de facto 
bound to a small range around zero. 

Explicit expressions for the potentials are rather involved [5] and will be 
given below in (16), considering again symmetric configurations only for the 
sake of simplicity (see [5, 6] for the general case). Note in there r and R are 
geometrical parameters (and m, I represent masses and inertia moments), while 
the dynamical variables are {ipn(t),$ n (t)}, n€ Z. 

Before writing these explicit expressions, we state that Ut and U s correspond 
to the torsional energy of the backbone and the base stacking energy respec- 
tively, and have simple harmonic expressions 11 in terms of the three-dimensional 
coordinates of involved nucleotide elements (albeit a rather involved expression 
in terms of torsion angles); as for the pairing interaction U p , this corresponds 
again to a harmonic potential in three-dimensional distance between bases in 
a pair. This latter choice, albeit not natural physically (a Morse potential as 
in the PB or BCP models would perhaps be more appropriate; see [24] for a 
discussion of different intrapair potentials within the Y model) was made in 
order to ease comparison with results obtained for the simple Y model, i.e. in 
order to focus attention on the improvements which have to be ascribed to the 
geometry of the model. Here are the explicit expressions for T and the Uc 



Ut = K t £„ [i - cos (K+i -#n)\ ; 

U s = K s J2n 2 [R 2 + r2 + 

- i? 2 COs(#„+i - tin) - r 2 COs[(z9„+i - l?„) + (fn+1 - <fin)] + 

- Rr (cos[(i? n+ i - ■!?„) + ip n+ i] + cos[(i9„+i - i? n ) - ip n ]) + 
+ Rr (cos(p n+ i) + cos(<p n ))] ; 

U p = (1/2) K p J2 n (a n - pf , with 

a 2 = 4 [a 2 + R 2 + d\ + 2 (Rd h cos ip n - aR cos & n - ad h cos((^ n + ??„))] . 



Here m is the base mass, / the momentum of inertia of the disk modelling the 
backbone units and and K u K Sl K p are, respectively the torsional, stacking and 
pairing coupling constants; <r„ represents the distance between bases in a pair 
(between end points of double pendulums), so that U p is a harmonic potential in 
the physical distance, albeit expressed by a non-harmonic function when dealing 
with angular variables. 

11 The possible relevance of nonlinear stacking interactions has been recently noted in [44]. 



T 



3 E„ mr 2 (Lp n f + 2mr(r + Rcos((p n ))i3 n Lp n + 
+ (I + m(R 2 +r 2 ) + 2mRr cos(^„)) U n ) 2 ; 




(16) 
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We work in the so-called contact approximation in which the equilibrium 
distance between the two basis, po vanishes (see Figure (2). In this approxima- 
tion the parameter a appearing above is given by a = R + dh- Proceeding as in 
the PB, Y and BCP cases above, we obtain field equations for the interpolating 
fields t) and Q(x,t); these are rather lengthy and reported explicitly in 
Appendix B of [5]. As we are mainly interested in soliton equations, we will 
introduce the travelling wave ansatz 

®(x,t) = $(x - vt) = <j>(z) , Q(x,t) = Q(x-vt) = 6(z) . (17) 

The equations for travelling wave (TW) solutions are 

pr 2 (j)" + fir(r + R cos 4>)0" = 

= -2aK p (a - R) sin (<£ + 0) + K s 5 2 Rr sin (0) {B'f + 
-R sm(<f>) {-2K p (a - R) + mrv 2 (9') 2 ) - ^ ; (18) 
pr(r + R cos </>) <f>" + [J + p(R 2 + r 2 + 2Rr cos <j>) 9" = 

= -2aK p (Rsm6 + {a- R)sm((l) + 0)) + (iRr sm(4>)[(4>') 2 + 2<j)'6'] . 

Here we have simplified the notation by introducing the constants 

/i = (mi) 2 - K s 5 2 ) , J = (Iv 2 - K t 5 2 ) . (19) 

It should be stressed that the ODEs (18) are obtained as a reduction of 
the nonlinear wave PDEs for $ and 6; the latter should be supplemented by 
boundary conditions. Requiring that the solutions go to an equilibrium for 
x — ► ±oo (and thus have finite energy), these are <fr(±oo,i) = <& x (±oo,t) = 
$t(±oo,t) = 0, 6(±oo,i) = 2n±7r, 6 x (±oo,i) = 6 t (±oo,i) = 0. These entail 
side conditions for (j>(z) and 0(z), i.e. 

^(±oo) = , 6»(±oo) = 2n±7r ; 0'(±oo) = , 0'(±oo) = . (20) 

We also note that the constraining potential U c makes that the fields (O, $), 
which in principles take values in S 1 x S 1 , are actually taking values in S* 1 x I 
(where I = (—A, A) is a real interval centered in zero, with A -C 7r); thus O is 
a topological field, while $ is a non-topological one (that is why the boundary 
conditions for $ do not allow nontrivial multiples of 2ir). 

2.2 Physical values of the parameters 

One of the nice and striking features of our composite model is that it supports 
solitonic solutions within a fully realistic range of all the physical parameters 
characterizing the DNA; this should be compared with the situation for sim- 
ple torsional models mentioned above, where unphysical coupling constants are 
needed to fit some experimental data. 

Let us therefore briefly discuss how the values of the parameters appearing 
in Eq. (16) are fixed. 
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R 


r 


dh 


a 


3.1 A 


2.7 A 


4.4 A 


7.5 A 



Table 1: The geometrical parameters in the composite DNA model. 





K s 




130 KJ/mol 


16.6 N/m 


3.5 N/m 



Table 2: The dynamical parameters in the composite DNA model. 



There are basically two types of parameters: kinematical ( the geometri- 
cal parameters R 1 r 1 dh,a 1 the mass m and the momentum of inertia I), and 
dynamical (the elastic coupling constants K tl K Sl K p ) 12 

The kinematical parameters can be evaluated by considering the chemical 
structure and the geometry of the DNA molecule. The order of magnitude of 
the mean value of the bases mass is m ~ 130 (atomic units). For the momentum 
of inertia of the disk we have I ~ 5 x 10~ 44 Kgm 2 , whereas the values of the 
geometrical parameters are given in Table 1. 

The determination of the numerical values of the coupling constants charac- 
terizing our model is more involved. Their order of magnitude can be estimated 
by considering the typical energy of hydrogen bonds (K p ) and the experimental 
results for the torsional rigidity of the DNA chain (K s and K t ) [5]. The results 
are given in Table 2. 

These values allow in particular to estimate the speed of torsion waves as- 
sociated to base torsion; this will be the speed of transverse elastic waves along 
the double chain. 

With the above values, the speed of elastic waves is estimated to be v s = 
Sy^Kg/m w 3Km/s; this is of the right order of magnitude, and well compatible 
with experimental data [32, 59]. 

Let us stress that the geometry of our composite model makes that using 
natural parameters one obtains predictions that nicely fits with the estimates of 
the structural properties and binding energies of the DNA: the induced optical 
frequencies and phonon speeds are of the same order of magnitude of those 
experimentally observed. This does not happen in simpler models, as stressed 
above when discussing the simple Y model. 

2.3 Nonlinear dynamics and soliton solutions 

Equations (18) are a system of two coupled, non linear ODEs. In general it 
cannot be solved analytically in closed form. One has to resort to numerical 

1 Notice that the values of the physical parameters given in this note differ slightly from 
those of Rcf ( [5] ) . The new values improve the estimates of the parameters but do not change 
the qualitative behavior of the model. 
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calculations in order to show that the system admits solutions satisfying the 
boundary conditions (20) [5] (see also Ref. ([14]). 

We note that the relevance of the dynamical system (18) goes well beyond 
DNA torsional dynamics: the same kind of equations appear in more general 
cases. 

In fact, we can as well consider (18) as describing the continuum limit of the 
torsional dynamics of a single molecular chain made of a disk and a pendulum. 
In this case, the pairing interaction for the DNA double chain is replaced by an 
external potential: 

V= -Ar 2 K p (cos 9 + cos((p + 9) - i cos y> - 3/2^ , (21) 

whereas the stacking and torsional interaction generate the ^-derivative terms. 

To further simplify our model - and the resulting explicit formulas - and 
concentrate on essential features of interest also beyond DNA dynamics, we also 
set R = r in (18). The resulting equations take the much simpler form 

H4>" + fj,{l + cos<}>)9" = 

= -AK p sin (<f> + 6) - n sin (<j>) {9'f + 2K p sm(</>) - ; 
M(l + cos </>)</>" + [(J/r 2 ) + 2/i(l + cos0)] 6" = 

= -AKp (sin 9 + sin(0 + 9)) + fi sin(0)[(0') 2 + 2<j>'e'\ , 

The previous form of the equations of motion will be taken as starting point 
to discuss two general features of the dynamical system. As we will see in 
the next two sections these features (the existence of a mechanism to select the 
speed of solitonic solutions and the slaving of the field <p) represent quite general 
consequences of nonlinear dynamics. We expect they will have a quite broad 
field of application in the context of non-linear Physics and Mechanics. 

Although in the general case one can find a solitonic solution of the system 
(22) only numerically, there is a particular case which admits analytical solu- 
tions. This is obtained by freezing the angle cj>, i.e by setting cf) = 0; note that 
if we force </>(z) — 0, we are actually considering a chain of simple pendulums, 
i.e. a sine-Gordon equation. 

This constraint can be accommodated in our setting in a dynamical way, by 
acting on the confining potential U c : this should be made stronger and stronger 
and the maximum angle cf>o will become smaller and smaller. 

Setting = and using (dU c /d<p)(0) = 0, the system (22) is equivalent to 
the equation 

fi9" = -2K p sm9 ; (23) 

the compatibility condition between the two equations of the system (22) is now 
given by 

J = 0. (24) 

Equation (23) has to be integrated with the boundary conditions (20). When 
\l < and n = 1, we have the kink 

6» = 4arctan[e /3z ], = 0, (25) 
13 
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Figure 3: The n = 1 soliton for the composite model of Figure 2. We plot in 
thick curves the fields 9 (on the left) and <j> (on the right) as a function of x 
for the static (v = 0) solution with physical values for the model parameters 
(reproduced from [6]). The solution is very similar to the corresponding one 
obtained for the Y model for what concerns the 9 field. 



where we have written, as in (12), (3 = y ~zjt- The solution (25) is of course 
the same as (11), i.e. the solution found in the context of the simple Y model. 

As for general solutions, we note these will still be indexed by the topological 
index n (which refers to the 9 behavior); the full equations (22) cannot be solved 
analytically in the general case (see below for a perturbativc approach), but they 
can be studied and solved numerically; the solution for n = 1 is displayed in 
Figure 3. 

2.4 Discussion 

It is quite clear that a weak point of this model is represented by the choice of 
the exceedingly simple potential U p ; and also by the the simplifying assumption 
p = (see (16) and comments thereafter). 

As mentioned above, this choice was justified by the will to ease comparison 
with results obtained with the simple Y model, i.e. to be able to focus on new 
features depending only on the more articulated geometry of the model. 13 

Needless to say, one should then consider the same model with more realistic 
pairing potentials - e.g. with the Morse potential used in the PB and BCP 
models (this is being done [14], and yields quite interesting preliminary results). 

It should be mentioned, in this respect, that investigations conducted within 
the framework of the simple Y model have shown that while dispersion relations 
are of course strongly affected by the contact approximation and by the choice of 
the pairing potential, these have very little effect on the soliton equations (pro- 
vided parameters in the pairing potential are set obeying to the same physical 
argument and considerations); this is shown in Figure 4. 

13 And also - in the present case - to be able to discuss interesting phenomena without being 
forced to tackle technically hard computations which could hide the physical and mechanical 
meaning of the results. 
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Figure 4: Comparison of the standard Yakushevich n = 1 soliton (11) (dotted 
curves) with the n = 1 solitons (solid curves) obtained numerically by relaxing 
the £o = contact approximation (left, from [23]) and with a Morse rather 
than harmonic intrapair potential (right, from [24]). This clearly shows solitons 
solutions are very little affected by the simplifying choices used in the Y model. 

3 Soliton's speed selection 

A general feature of sine-Gordon solitons (and more in general of relativistic 
solitons) is that the soliton speed \v\ is a free parameter, which can be fixed 
by choosing initial conditions and is bounded from above by a limiting value 
co. This is a consequence of the Lorentz symmetry of the equation that fixes 
a limiting upper bound for the speed, whereas v can be changed by applying a 
boost. 

Very often it happens that non-linear systems (e.g the DNA chain, but also 
reaction-diffusion equations [38], tsunami equations [34], etc.) that allow for 
solitary, non dispersive excitations, somehow select the order of magnitude for 
the speed of propagation of this excitation, in the sense that when solitons are 
experimentally observed, they turn out to have a speed of a well defined order 
of magnitude. 14 

Thus, in practical situations, whenever the experiments give a well-defined 
value for the propagation speed of the soliton, the speed degeneracy represents 
a loss of predictive power of the model. It is quite remarkable that our model 
has a built-in mechanism for selecting the soliton speed; it is essential for this 
mechanism that we have (at least) a two-components system. 

3.1 Speed selection in the composite DNA model 

One can easily realize that the compatibility condition (24) fixes the speed of 
propagation of the soliton (25) to the speed c t of the transverse sound waves 

14 Wc observe that if soliton excitations are relevant in DNA transcription, they should also 
have some built-in speed selection mechanism: in fact, the order of magnitude of the speed of 
the transcription bubble along the DNA double helix is well defined (and such to coordinate 
with the synthesis of RNA messenger by RNA Polymerase), as mentioned in the caption to 
Figure 1. 
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supported by the elastic torsional forces acting on the disk 



v = ct = V^/I, (26) 

where to t — K t S 2 . 

Moreover, as the soliton exist only for \i < 0, the soliton speed is bounded 
from above by the speed c s of transverse sound waves supported by the elastic 
stacking forces acting on the pendulum, i.e. 



v < c s = ^Jujs/m, (27) 

where w s = K s 6 2 . In view of Eq. (26), this implies a constraint on the stacking, 
torsional coupling constants and kinematical parameters of the system: 

K t K s 

-r < — . 28 
I m 

This selection mechanism for the soliton speed gives a nice and simple way to 
produce solitons with a given speed in double pendulums molecular chains. To 
select the soliton speed one just needs to tune the torsional and stacking coupling 
constants and the kinematical parameters of the chain such that Eqs. (26) and 
(28) are satisfied. Acting on the confining potential U c , making it stronger and 
stronger, one obtains the single pendulum limit of the double pendulums chain. 
The angle is frozen to <j> = and a SG soliton with a speed equal to that of the 
transverse sound waves supported by the torsional forces acting on the disk is 
selected. 15 

Notice that the mechanism described here can be obviously used to devise 
and realize non-linear media where solitons propagate at a given fixed speed. 



3.2 Speed selection in general models 

The mechanism for selecting the soliton speed described above for the molecular 
chain model (22) is rather generic. It is related to the existence of a condition- 
ally conserved quantity J 16 and it is rather independent from the specific form 
of the interactions characterizing the model. The proposed mechanism will 
work whenever we have a nonlinear mechanical system satisfying some general 
conditions: 

1. The system must have at least two degrees of freedom {X, Y) at each site, 
which in the continuum approximation will give two interpolating fields 
X(x,t), Y(x,t) and are characterized by masses (or moments of inertia) 
to, M with to ^ M; 

15 Thus in a way this mechanism is related to the fact the simple pendulum model is not 
structurally stable, and should be seen as the singular limit of a class of more general systems. 
Note this class could be not unique: e.g. for the pendulum case, one could consider a chain 
of coupled pendulums made of elastic beams, obtaining the sine-Gordon equation as the limit 
case when beams are infinitely rigid. 

16 This is the momentum conjugate to the angle 9 evaluated at <j> = [7], J = (dC/dO)<f, = o. 
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2. There should be at least two types of interactions: (a) An elastic force 
(coupling constant K t ) originated by the interaction between neighboring 
sites on the chain; (6) A non linear force (coupling constant K p ) acting on 
the single site; 

3. There should be a confining potential U c that limits the range of variation 
of one degree of freedom (e.g V) and allows to freeze Y. That is, by making 
U c steeper and steeper we implement the dynamical reduction from two 
degrees of freedom to one single degree of freedom X; 

4. Freezing the degree of freedom Y we have both a conservation law for the 
momenta conjugate to X and solitonic solutions for X. 

4 Perturbative expansion and slaving 

A drawback of the general composite Y models (16) is that the equations de- 
scribing its dynamics are far too complex to be exactly solved. 

On the other hand, the doubling of degrees of freedom introduces a natural 
separation between topological and non-topological degrees of freedom. This 
separation, together with the fact that the composite model allows for the exact 
Y soliton (25) when the non-topological field <fi is frozen, opens the way to a 
perturbative treatment of our model. 

The simplest way to deal with our non linear system at the perturbative 
level, is to consider it as embedded in a family of double pendulums chains, 
which has a single pendulum chain as a special case. 

We will first consider a chain of simple pendulums (of length R and mass M); 
we will then look for solutions of the double pendulums model by perturbing 
the system near the single pendulum Y solutions (25). 

As the simple pendulum limit of our double pendulums chain involves a 
reduction of the number of degree of freedom, we will have to deal with a 
singular perturbative expansion. 

A way to obtain the simple pendulums chain (and the sine-Gordon equation 
in the continuum limit) from our model is to let the length of second pendulums 
go to zero. In this case our set of parameters becomes redundant, as the positions 
of the masses of the pendulums coincide in space, so that only the total mass is 
relevant. A similar argument can be used to show that only the total coupling 
strength K = K s + K t is relevant. 

If the double pendulums chain is seen as a (singular) perturbation of the 
simple pendulums, one is naturally led to look for travelling wave solutions as 
perturbations of the standard sine-Gordon solitons. 

This means we look for solutions to the equations (18) in the form of a series 
expansion in a small parameter e. We will correspondingly also expand in the 
same parameter the parameters appearing in the model and in the solution: the 
geometrical parameters, the masses appearing in our model, the two coupling 
constants K t and K s and also allow for modification of the speed by expanding 
it as well [30, 57]. 
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Therefore, the series expansion we adopt are as follows: 



m 



v 



M 



R 



K t 
K s 



d 



r 



i?o + e{89) 

e{54>) 

e(5r) 

A - e{Sr) 

e{5m) 

M tot - e(5m) 

v a + e(5v) 

e(6K) 

K - e(5K) 



i?o + e^i +s 2 -d 2 + .... , 

e(pi + e 2 4>2 + .... ; 

eri + e 2 r 2 + ... , 

A — er\ — e 2 r 2 — ... ; 

smi + e 2 m 2 + ... , 

M tot - m - smi - e 2 m 2 - ... ; 

v + svi + e 2 v 2 + ... ; 

efci + e 2 k 2 + ... , 

K — eki — e 2 k 2 — ... . 



(29) 



where $0 and 4>q are the limiting simple pendulum solutions given by (25). 

Inserting the series expansions (29) in the equations for TW solutions (18) 
we will obtain perturbative solutions of our non linear dynamical system at the 
various order of the perturbative expansion in the parameter e. 

We will skip here the computational details, which can be found in Ref. 
[6, 8]; the results one obtains in this way are summarized as follows. 

• At zero order of the perturbation theory the single pendulum solitonic 
solution (25) is reproduced. 

• First and higher order corrections to i9o and <fio can be explicitly calcu- 
lated. They exhibit the following striking feature: the non-topological 
field (j)k turns out to be completely determined algebraically - we say 
then it is slaved - order by order by the topological one 9k- Thus flk 
is obtained as the solution to a differential equation which depends on 
{0o, <Afe-i; ^0 • • • #fc-i}, while (f>k is determined algebraically (no differ- 
ential equation involved!) by {0o, <Afc-i; $0 ■ ■ ■ #fc-i}- 

• One can expand £ in a series in e, C = J2k £k ^k' : then at any order in the 
perturbative expansion the Lagrangian Ck depends on 4>k but it is inde- 
pendent of the momentum conjugated to 4>k- Thus 4>k can be considered 
as an auxiliary field, entering in the Lagrangian only algebraically (and 
not differentially). We can express this fact by saying that the field is a 
auxiliary field in perturbation. 

Once again these remarkable features are not specific to the DNA model con- 
sidered here, but do quite obviously apply to a much wider class of mechanical 
(and field-theoretical) models. 

Actually, one can state that whenever a two-component evolutionary equa- 
tion can be expanded in series so that one of the two fields is an auxiliary field in 
perturbation, then it will be slaved in the perturbative expansion, and pertur- 
bative solutions will admit the solution of the simpler PDE obtained by freezing 
the field which is auxiliary in perturbation to zero as the limit for e — * 0. 
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5 Solitons in polyethylene crystals 



In this section we will substantiate our claim that the mechanisms (in partic- 
ular, the soliton speed selection and the slaving) devised in the study of DNA 
are significant for more general classes of polymeric macromolecules by show- 
ing how they apply to a Polyethylene (PE in the following) chain in a crystal 
environment. 

The (nonlinear) dynamics of crystalline PE chains has been discussed in 
various papers [3, 29, 36, 46, 47, 60, 61, 62] following pioneering modelling 
work by Kirkwood [33] and by Mansfield and Boyd [37]; here we will follow the 
approach by Zhang and Collins [62]. 

5.1 General features 

It has been argued that the twisting and the elongation (or compression) of the 
PE chain can be described by elementary excitations called twistons, smooth 
twists of the PE chain accompanied with a contraction (or elongation) of the 
CH 2 units of the chain [62]. These twiston excitations may be relevant to de- 
scribe the propagation of conformational defects in the PE chain or, more in 
general, its molecular dynamics. Using realistic intramolecular and intermolecu- 
lar interactions, it has been shown that twisting and elongation of the PE chain 
can be described by several types of topological sine-Gordon solitons [62]. This 
is essentially due to the symmetry of the crystals structure, which allows for 
many equivalent ground states of the interchain potential. As in the continuum 
limit the dynamics of the PE chain is modelled by coupled equations of (mul- 
tiple) sine-Gordon type, we expect that the speed selection discussed in section 
3 and the slaving discussed in section 4 apply also to the case of crystalline 
PE chain. Before discussing these effects let us briefly summarize the relevant 
features of the Zhang-Collins PE model [62]. 

Each CH 2 unit of the PE chain is considered as a rigid group with mass 
m; it is labeled by an index n and has both intermolecular interactions with 
the whole crystal environment, modelled by a potential U a , and intramolecular 
interactions with its next-neighbor CH2 units, modelled by a potential Up. 

Owing to the symmetry of the crystal environment a cylindrical coordinate 
system is appropriate. The position of the n-th carbon atom is given by the 
cylindrical coordinates {p n , <fn,Cn)- 

The dynamics of the PE chain is described by the Lagrangian 

£ = T - U a - Up , (30) 
with kinetic energy given by 

n 

The intermolecular potential energy U a is the sum of the effective potential U 
for each CH 2 unit [62]: U a — J2 n U{<Pn, Cn)- Owing to the symmetry of the 
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crystal, U has a 7r-periodic part Uo(ip n ) and a 27r-periodic part Ui((p n ,(n), 

U(y n ,(n) = U (tp n ) + B sin(<^„ - <p n ) cos(ttC„/c) (32) 
where B, <p, c are constants and, neglecting a small constant phase, we have 

Uq = — A cos (2 (p n ) — A cos(4(^„). (33) 

Note that U(ip, £) in the whole is thus 27r-periodic in ip and 7r£/c. 

The intramolecular interaction Up can be given in terms of simple sums of 
the energy of C-C bonds, C-C-C bond bends and C-C-C-C torsions [62]. 

In order to describe the soliton's speed selection and the slaving in the PE 
chain we will consider a simplified model for the PE dynamics. Note that since 
the variation of p n is rather small even in the twisting region [62] , one can treat 
p n as a constant, at least in a first approximation. We will therefore set p n = r U; 
identically, in the Lagrangian (30). This assumption allows us to get rid of one 
degree of freedom, which is not essential for our purposes. 

We will also neglect the term in the potential Uq which is 7r/2-periodic; 
that is, we set A = in Eq. (33). This approximation allows us to reduce our 
Lagrangian to that describing a set of coupled double sine-Gordon equations. As 
well known, the double sine-Gordon equation is also intcgrable [10] and admits 
topological soliton solutions, so that it has to be expected that the general 
mechanisms at work in the simple sine-Gordon case would also apply here. 17 

Note that for A = 0, the potential (32) has a minimum for £„ — and (p n 
identified by the condition 

B cos{ip n - ip n ) cos(7rC„/c) + 2 A sin(2<^„) = 0. (34) 

In order to keep the calculations as simple as possible we will determine the 
phase ip n so that U has a minimum for = 0, ip n = —ir/2 (in terms of the 
fields introduced in (36) below, this corresponds to "J = $ = 0). This 

condition provides <p = 0. 18 

In the continuum limit n — > x/c, where x gives now the coordinate along the 
chain axis, the system is described by the Lagrangian density 

Tfi ' 7TU 

C = — (r 2 ip 2 +ufj - -(fc 2 c 2 ^ + k 3 c 2 u 2 x ) - EsinV>cos( — ) + ,4 cos 2-0, (35) 



17 In this sense, a rougher approximation reducing the problem to coupled (simple) sine- 
Gordon equations would also not change the qualitative features of the problem, while al- 
lowing a considerable simplification of its analytical treatment. As shown below, albeit this 
statement is correct in general, with the physical values of the PE model we are considering, 
a certain compatibility condition applying in the sine-Gordon approximation is not satisfied 
(and solitons like those we arc interested in cannot exist), while the equivalent condition is 
satisfied in the double sine-Gordon model. 

18 The value of <p n given in Ref. [62], is slightly different, <p n = 0.287T. Note however that 
this refers to the general model with 4^0, and phases are chosen in [62] so to have minima 
of the potential for ( n = ip n = 0; thus we are applying the same criterion used in [62] for the 
general model to our simplified case A = 0. 



20 



where the fields ip(x, t), u(x, t) represent the continuum limit of the displacement 
of the coordinates ip n {t), Cn(t) from their equilibrium values, tq is the equilibrium 
value of p, and the constants k 2 , k 3 depend both on the elastic coupling constants 
of the intramolecular potential and on the geometrical parameters of the PE 
chain [62]. 

Using suitable approximations, one can derive twiston solutions for the 
Euler-Lagrange field equations stemming from the Lagrangian (35). These so- 
lutions have the form of a kink for the field ip, topological sine-Gordon soliton 
for the field ttu/c [62]. These solutions bear a strong resemblance with the so- 
lutions of the composite DNA model described in the previous section. This 
is not surprising, as the models themselves arc very similar; in particular, the 
Zhang-Collins model for polyethylene would translate into a "composite" model 
of Cocco-Monasson type [11] model for DNA. 

It is therefore natural to look for a soliton speed selection mechanism, and 
a slaving mechanism, similar to those discussed in section 3 and in section 4. 

5.2 Soliton's speed selection in the PE chain 
5.2.1 The double-sine-Gordon model 

Performing in Eq. (35) the field redefinitions 

^ = tf-£, u=-$, (36) 

the Lagrangian becomes 

7TI C 2 1 

£ = !7(r^ t 2 + ^$ t 2 )--(fc2C 2 ^ + fc3— $2) + Bcos ^ cos $_ j4cos2 ^ . (37) 

Notice that now both *S? and $ can be considered as angular coordinates with < 
$ < 2n. Introducing the symmetric and antisymmetric field combinations 

X+ = (* + *), X- = (*-*), (38) 

the field equations for travelling waves solutions x+(£)> X-(0 (where £ 
give 

A 1 X+" + JX-" + 2Btt 2 sinx+-4A^ 2 sin(x++X-) = 0, 
A 4 X-" + JX+" + 2Stt 2 sinx--4yL7r 2 sin(x++X-) = 0, 

where the prime denotes derivation with respect to £ and p, J are given in terms 
of the elastic and geometric parameters of the model by 

p = mv 2 (-K 2 rl +c 2 ) - c 2 (k 2 n 2 + k 3 c 2 ) , , . 

J = mv 2 (TT 2 r 2 - c 2 ) - c 2 (k 2 n 2 - k 3 c 2 ) . [W) 

The system (39) gives a nice and simple description of the propagation of 
twistons in terms of symmetric and antisymmetric combinations of twisting and 
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m 


= 14.1 g/mol 


fc 2 


= 109.263 kJ/mol 


ra 


= 0.4236 A 


h 


= 1860.655 kJ/(A 2 mol) 


c 


= 1.274 A 


B 


= 1.52 kJ/mol 


A 


= 0.32 kJ/mol 




= 6.33 



Table 3: The dynamical, geometrical and coupling parameters in the Zhang- 
Collins polyethylene model, as given in [62]. 



elongation modes in the PE chain. It has the form of two coupled equations of 
sine- Gordon type. 

The system has many degenerate ground states \- = 2niir, x+ = 2n 2 7r, 
with ni,ri2 integers. We therefore expect the existence of topological sine- 
Gordon solitons connecting all these vacua. We will not discuss here the general 
solutions of the system (39), but just point out that the system gives a very 
simple realization of the speed selection mechanism discussed in Sect. 3. 

If wc are interested in the symmetric solutions of the system we can set 
in Eqs. (39) \- — 0- This can be realized dynamically introducing in the 
Lagrangian (37) a confining potential U c (x-) satisfying (dU c /dx-)o = and 
rising (sharply) with |x_|, so to freeze the degree of freedom \- to its vacuum 
configuration \- — (i.e. it forces symmetric configurations, <I> = <3>). It is 
easily seen that for \- = the field equations (39) become 

MX +" +( 2B-4A)ir 2 sin X+ - 0, (41) 

whereas the compatibility condition between the two equations (39) reads 

JB + 2A(fi-J) = . (42) 

Comparing Eq. (41) with Eq. (23) one easily realizes that for n/(B — 2A) < 
0, Eq. (41) admits SG solitonic solutions given by (25). The compatibility 
condition (42) fixes the speed of the soliton 



/ &27T 2 W - fc 3 C 2 

V ( s +) = ±cJ 2 2-, 43 



where a; is a dimcnsionless parameter given by 

For the values of parameters A, B proposed by Zhang and Collins [62] and 
reported in Table 3 we have u> = 6.33. Moreover, from these values of the 
parameters A, B it follows B > 2A, so that the condition for the existence of 
solitonic solutions becomes [i < 0. 

This condition implies a maximum speed for the soliton 



/ fee 2 + k 2 Tr 2 

< V M = ca — 2-2T • 45 ) 
V m(c 2 + 7r 2 rg) 
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Needless to say, forcing \+ — we would obtain a fixing of the soliton solu- 
tion speed for antisymmetric fields provided by the same compatibility condition 

(42) , so that Eqs. (43) and (45) still hold for the antisymmetric solution. 
Note that the argument of the square root in (43) could be negative, de- 
pending on the values of the parameters. When this happens, no symmetric or 
antisymmetric travelling wave solution is possible. Moreover, the soliton speed 

(43) has to be smaller than the maximum allowed speed (45). 

Let us now check the existence of the symmetric and antisymmetric soliton 
solutions using realistic values for the physical parameters of the model. The 
values of the physical parameters characterizing the PE chain [62] arc collected 
in Table 3; they indicate that kiT^u is much bigger (at least two times) than 
fc 3 c 2 , whereas 7rr w > c, allowing for both symmetric and antisymmetric solitons 
with fixed speed 

v (±) - 7 x 10 3 m/s. (46) 

On the other hand the maximum soliton speed (45) turns out to be well 
above the speed of the symmetric or antisymmetric solutions. From the values 
of the parameters of Table 3 we get 

V M ^ 1-2 x 10 4 m/s. (47) 
5.2.2 The sine-Gordon approximation 

We would like to stress that the soliton's speed selection mechanism is a general 
and rather robust effect, which does not depend on the details of the model we 
are considering - but the very existence of soliton may depend on the value of 
physical parameters. To have a flavor of this fact let us consider (as anticipated 
in a footnote above) an even more simplified model for the PE. 

We will neglect completely the 7r-periodic part in the potential U (tp, u) in 
(32), i.e. we will set there U a = 0; note this does not change the overall period- 
icity properties of U a . This is a very rough approximation, but it allows us to 
reduce our Lagrangian (35) to that describing a simple two-fields sine-Gordon 
model. This situation represents just the particular case A = {u> = 1) of our 
general equations. 

By forcing the degree of freedom X- to its vacuum configuration X- = 
(the same discussion would apply interchanging the roles of x+ and X-)i the 
field equations (41) become now 

MX+" + 2Btt 2 sinx+ = , (48) 

whereas the compatibility condition (42) simplifies to J = 0. This selects the 
soliton speed, 

v = Vfo = ± M fe^. (49) 
w V m(c 2 - 7r 2 rg) 

The values of the physical parameters of Table 3 indicate that k$c is much 
bigger (at least three times) than k^ 2 . Equation (49) requires therefore c > nr . 
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The condition for the existence of solitonic solutions (45) together with the 
previous equation imply that the soliton solutions of our simplified model can 
exist only for r\ < (^2/^3). The values of the parameters reported in Table 3 
do not satisfy this condition. Thus for the values of the parameters proposed 
by Zhang and Collins, the symmetric and antisymmetric soliton solutions (41) 
of the more general model (35) are allowed, but the rougher approximation 
obtained setting A = is not allowing such solutions. 



5.3 Speed selection and conditionally conserved quantities 

In Ref. [7] it has been pointed out that the speed selection mechanism is related 
to the existence of a conditionally conserved quantity. It is not difficult to 
identify this quantity for the PE model we are discussing here. For the case 
in which a confining potential for the field X- i s introduced, from the field 
equations (39) one can easily derive the following equation 

F + , (50) 



where 



W+ = (BJ + 2A(fi - J)) X ' + + + 2A( M - J)) X '_, 
F+= -2^B[(2A-B)sin(x-) + 4AB(sin(x++X-)-sin X+ )] (51) 
-[dU c ( X -)/d X -]. 

If X- = 0, then F + ( X - = 0) = and W+ is conserved. Moreover, when x- = 0, 
the conserved quantity becomes W + = (BJ+2A(p,— J))x'+, which is compatible 
with the existence of solitonic solutions to (41) only for BJ + 2A([i — J) = 0, 
i.e. only when the soliton speed is given by (43). 

Similarly, when a confining potential U c (x+) for the field x+ is introduced, 
we have the conditionally conserved quantity W- and a force F_, whose ex- 
pressions are obtained by interchanging x+ an d X- m (51)- Again, when x+ 
is frozen to its vacuum value x+ = 0, W- is conserved and the speed of the 
X--soliton (if this is allowed) is fixed to the value (49). 



5.4 Series expansion, and slaving 

In order to discuss series expansion around a given solution, we will use the 
coordinates x± , and for notational convenience we rewrite these as 

P = x+ , Q = x- ; (52) 

moreover, we will work directly in the space of functions of the variable £ = 
(x±vt). 
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5.4.1 The double-sine-Gordon model 



In this notation, and denoting again £ derivatives by a prime, the Lagrangian 
(37) can be written as 



M(P' 2 + Q' 2 + 2SP'Q') +2B(cosP+cosQ) - cos(P+Q) , (53) 



where we have further simplified the notation by writing 

M=4, 5=4- (54) 

We will consider a Lagrangian £ £ depending on a parameter e so that for 
e — > it reduces to a (sine-Gordon) Lagrangian Co which only depends on 
(P, P') plus a term constraining Q = 0, while for e = 1 it is just the Lagrangian 
(53). The simplest Lagrangian with these properties is 

C e = i M(P' 2 + sQ' 2 ) + 2SeP'Q' + 2B(cosP + cosQ) - AAecos{P + Q) . 

(55) 

The corresponding Euler-Lagrange equations are 

MP" + eSQ" + 2Bsm(P) + 4 As sin(P + Q) = 0, . . 

eMQ" + eSP" + 25sin(Q) + 4Assin(P + Q) = 0. ^ ' 

We will then consider e series expansions both for the functions of £ and for 
the parameters; thus we write 

P= P +ePi+e 2 P 2 +e 3 P 3 + ... , 

Q= Q +e 2 (?2 + e 3 (?3 + •■■ ; 

A= A +eAi +e 2 ^ 2 +e 3 A 3 + ... , . . 

B= B + sB 1 +e 2 B 2 + e 3 B 3 + ... ; 

M = M + eMi + £ 2 M 2 + e 3 M 3 + ... , 

,5= ,So + eS!+ e 2 ,S 2 + e 3 5 3 + ... . 

Note we are also expanding in series the constants B, M, S, as prescribed by the 
general Poincare-Lindestedt procedure. 19 
We now expand C £ in e, 

C e = L + eL 1 +e 2 L 2 +e 3 L 3 + ... . (58) 



19 Moreover, M and S depend not only on the couplings and the geometric parameters of the 
model, but also on the wave speed v, which should in any case be allowed to vary with e. It 
may be worth remarking that the terms of the parameter expansions would be determined by 
imposing the projection of higher order terms in the expansion for P on the space of solutions 
to the sine-Gordon equation to vanish. As we arc not going to discuss this aspect (that is, 
implicitly, how the speed of the solution depends on e), but only in displaying slaving of the 
$ field — i.e. of the Q variable - these terms will remain undetermined here. The reader can 
easily check that slaving shows up as well by just setting (A, B, M, S) = (Aq, Bo, Mo, So). 
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At order zero, we get (as required) 

L = \ MoP'o + 2B cos(Po) + 2P cos(Q ), (59) 

which yields the sine-Gordon equation for P , 

P"o = - (2B /M Q ) sin(Po) , (60) 

together with sin(Qo) = 0; this in turn enforces (up to a suitable choice of the 
origin for $) 

Qo = . (61) 

At order e, we get 

Li = \ MiP'l + MoP'oP'i+2B 1 [l + cos(Po)}-2B Pism(P ) + 4A cos(P ) ; 

(62) 

the Euler-Lagrange equation with respect to Pi is identically satisfied (it just 
yields back the sine-Gordon equation for Po). 
At order e 2 , we get 

L 2 = (l/2)[M P'l + M 2 P' 2 } + MrP'oP'i + P' (M P' 2 + SoQ\) + 

- BoQj + (2B 2 - B P? + 4Ai) cos(P ) + 

- 2[BiP! + B P 2 + 2A Q (P 1 + Qx)] sin(Po) • 

(63) 

The equation issued by variation with respect to Pi, taking into account 
that Po solves (60), reads 

M P"i = - 2[BoM cos(Po)] Pi + {2/M Q )[{B Q Mi-BiM Q )-2AoM Q } sin(P ) • 

(64) 

As for the variation with respect to Qi (obviously taking again (60) into 
account), it yields 

Qi = (^-2^)sin ( P ). (65) 

This is just a relation between Qi and Po, and uniquely determines Qi in terms 
of the latter and of the parameters appearing in the equation. In particular, 
it is not a differential equation, and shows that Q (that is, the field <f> in the 
notation used earlier on) is slaved to P (that is, to the field *S?) at leading order. 

Actually, it is easy to see that the same will happen at all orders and not 
just at leading order. That is, the terms Pk will be determined by linear non- 
autonomous differential equations, the non-autonomous terms of these depend- 
ing on terms {Po, Pi, Pfe-i; Qi, Qk-i} of lower degree; on the other hand, 
the terms Qk will be determined algebraically in terms of the {Po, Pk}- 

Needless to say, the explicit equations (and the partial Lagrangians Lk too) 
will become quickly quite involved. Thus we will just give the equations satisfied 

20 Using this fact, the equations for the Pj. can be written in terms of the Pj (with j < k) 
alone, i.e. without explicit use of the Qj . 
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by e 2 order terms for P and Q. The term P 2 is obtained as solution to the 
differential equation 



P M 3 P" 2 = [(B MoS 2 - 2A M 2 S ) sin(P )] P'l + [B M cos(P )] P2 + 
- 2[B Mo{2A Mo + PiM - B M{) cos(P )] Pi + 
+ [P 2 M 2 sin(P )]P 2 + 2[(2A + B 1 )B M M 1 - {2A 1 + B 2 )B M 2 + 
+(M M 2 - M 2 )Bl + (B0S0 - 2A a M ) 2 cos(P )] sin(P ) ; 

(66) 

as for Q 2 , this is given by 



Q2 



(2B 2 M 2 ) 



__„ , LV 2A Mo(2(P 1 Mo + PoS )-Af 2 P' 2 ) + 

+ B (2B (M S 1 - M 1 S ) - ^M 2 + M 2 S P'l)) sin(P ) 
+ 2M (2A M - PoSo) cos(P )((2A - B ) sin(P ) - B Pi) 



(67) 



5.4.2 The sine-Gordon approximation 



Slaving of the Q degree of freedom is present also in the sine-Gordon approxi- 
mation considered above, i.e. for A = 0. In this case one obtains slightly simpler 
explicit formulas, easily obtained by setting A = A\ = in the previous ones. 

As recalled above, see the footnote following Eq.(57), if one is only interested 
in slaving a simpler series expansion - with A, B, M, S not depending on e - 
would also be possible. We will now adopt this simplified expansion with the 
sine-Gordon approximation (i.e. A = 0) in order to see the slaving mechanism 
at work for the PE model in its simplest setting. 

With the same notation as above, the Lagrangian (37) with A = reads 



M(P' 2 + Q' 2 + 2SP'Q') + 2P(cosP + cosQ) 



The Lagrangian C e can be chosen as 



1 



M(P' 2 + eQ' 2 ) + 2SeP'Q' + 2P(cos P + cos Q) 



and the corresponding Eulcr-Lagrange equations are 

MP" + eSQ" + 2Psin(P) = 0, 
sMQ" + eSP" + 2Psin(Q) = 0. 



(68) 



(69) 



(70) 



With the simplified series expansions, and expanding C e in e as before, we 
have partial Lagrangians Lk corresponding to terms of order e h . 

At order zero, we get of course (59) (except for Bo, Mo now reading B, M), 
which yields the sine-Gordon equation for Po, 



P"o - - (2P/M) sin(Po) , 



(71) 



together with Q = 0. 
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At order e, we get 



L x = MP' P\ - 2PP 1 sin(P ) , (72) 

with the Euler-Lagrange equation with respect to Pi identically satisfied (it just 
yields back the sine-Gordon equation for P ). 
At order e 2 , we have 

L 2 = (l/2)(MP'i) + P' (MP' 2 + SQ\) - BQ\ . (73) 

The equation issued by variation with respect to Pi yields, taking into account 
that P solves (71), 

MP" i = - 2(BMcos(P )) Pi . (74) 

As for the variation with respect to Qi, taking again (71) into account it 
yields 

M Qi = 2BS sin(Po) • (75) 

Thus again we get a relation between Qi and Po (not a differential equation), 
which uniquely determines Qi in terms of the latter and of the parameters 
appearing in the equation. 

Albeit simpler than in the general setting, explicit expressions still become 
quite involved at higher orders, and we just give the equations satisfied by e 2 
terms. As for P 2 , it solves the equation 

M 3 P" 2 = -2PM 2 cos(P )P 2 + [PM 2 P 2 + M5 2 P'o + 2BS 2 cos(P )] sin(P ) . 
The next-to-leading order term Q 2 for Q is given by 

Q 2 = (2PM 2 )- 1 [M 2 SP'l + 2PM5cos(P )(P 1 +sin(P ))] . (77) 

6 Conclusions and discussion 

Almost thirty years after the seminal paper by Englander, Kallenbach, Heeger, 
Krumhansl and Litwin, non linear mechanical models of DNA still represent 
an active area of research, and a toll for trying to tackle fundamental problems 
such as the denaturation and transcription processes. 

The nice feature of this mechanical approach, not shared by approaches using 
full molecular dynamics, lies in its simplicity. This simplicity allows to model 
general features of DNA and to extract relevant information with relatively few 
analytical and/or computational effort. 

Simple DNA mechanical models are obviously too simple to take into account 
the full complexity of the DNA macromolcculc; on the other hand, they may 
well be able to describe DNA dynamics for what pertains to specific biological 
processes - such as DNA thermal denaturation or the formation and dynamics of 
open bubbles to which RNA Polymerase could bind in the transcription process. 
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Simple DNA mechanical models as the ones formulated by Peyrard and 
Bishop and by Yakushevich were able to provide correct qualitative predictions, 
and fit the order of magnitude of biologically relevant and physically observable 
quantities (e.g. the frequency of small amplitude oscillations, characteristic 
scales of breathers and some statistical mechanics features [26, 39, 41, 42] in 
the denaturation transition, and the size of solitonic excitations [22, 26, 42, 54] 
in the context of DNA transcription) ; on the other hand they failed completely 
when confronted to other physical quantities which are directly observable in 
modern single-molecule experiments [35, 43] and which involve elastic properties 
of the DNA molecule, such as the speed of transverse elastic waves [55]. 

The new generation of DNA mechanical models, in which there are more 
than one degrees of freedom per nucleotide, such as the BCP and CM models 
in the context of DNA denaturation [1, 2, 11, 39, 40], and the composite model 
discussed in this note for what concerns torsional DNA dynamics and transcrip- 
tion, represent a big improvement in the direction of a more accurate modelling 
of DNA still retaining the attractive features of simple models. 

In fact, on the one hand they remain simple enough so that their dynamics 
can be at least controlled, if not completely solved, at analytical level; on the 
other hand they allow for a more realistic description of the DNA complexity. 

Focusing on torsional DNA dynamics and transcription, if one considers our 
composite model then the predicted speed for optical and sound excitation in 
the DNA chain fits the order of magnitude of the experimental data that cannot 
be fitted by the simple Y model with physically acceptable coupling constants. 

Moreover, the greater number of degrees of freedom per site - and more 
specifically the fact one of these refers to the homogeneous (sugar-phosphate 
backbone) of the DNA molecule, the other to the non-homogeneous (nitrogen 
Watson-Crick bases) of it - enables one to introduce in a natural way those 
inhomogencitics in the DNA chain (in the form of different basis sequences) 
that are necessary for the codification of the genetic information. 

We also stress that real DNA lives in a highly viscous (at the molecular scale) 
fluid and is subject to thermal noise. These features could be implemented more 
realistically within a composite model able to take into account the differences 
between the external (backbone) and the internal (bases) parts of the DNA 
molecule. 

Apart from the specific problems of DNA modelling, and maybe more rel- 
evantly for the general community working in nonlinear systems and nonlinear 
Mechanics, the research activity on DNA dynamical modelling has also con- 
tributed to focusing on previously unnoticed mechanisms and deepen our un- 
derstanding on nonlinear phenomena; it thus became a source of new ideas in 
the field. 21 

In a separate but related development, Saccomandi and Sgura [44] have 
realized that chains with fully nonlinear elastic nearest neighbor coupling would 
present peculiar features, in particular in these the solitonic excitations would 

21 For a discussion of this statement in relation to breather- type nonlinear excitations, we 
refer to the book of Peyrard and Dauxois [42]. 
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have a strictly finite size - thus be compactons rather than ordinary solitons (see 
also [49] in this respect). This mechanism can also be generalized and extended 
to more general systems than the DNA molecule [15, 25, 27]. 

For what concerns the model discussed here [5, 6], the mechanism for select- 
ing the speed of solitons by tuning the physical parameters of the system on 
the one hand [7], and the separation in slaving and master fields [8] described 
in this note are two nice examples of this statement. 

In particular, the speed selection mechanism which has been originally dis- 
covered for the DNA composite model [5, 7], can be generalized for an ample 
class of molecular chain models [7] and could find broad applications to de- 
vise and realize nonlinear media where solitary wave excitations propagate at a 
selected speed. 

Also the separation of the degrees of freedom in "master" and "slave" seems 
not to be limited to DNA non linear dynamics, but to be a quite generic feature 
of this ample class of nonlinear systems. It may be very useful for separating 
in a hierarchical way the different degrees of freedom that are relevant to the 
dynamics of the nonlinear system, and be a guiding principle to make easier 
perturbation analysis of such systems. 

We illustrated the possibility of a wider applications of these two mechanism 
in the study of (polymeric) macromolecules by considering a concrete different 
system, i.e. polyethylene, and showing by explicit computations how they apply 
in this case as well. 
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